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Abstract 

In  this  paper,  we  examine  the  effectiveness  of  absorbing  layers  as  non-reflecting  com¬ 
putational  boundaries  for  the  Euler  equations.  The  absorbing-layer  equations  are  simply 
obtained  by  splitting  the  governing  equations  in  the  coordinate  directions  and  introducing 
absorption  coefficients  in  each  split  equation.  This  methodology  is  similar  to  that  used  by 
Berenger  for  the  numerical  solutions  of  Maxwell’s  equations.  Specifically,  we  apply  this 
methodology  to  three  physical  problems  -  shock-vortex  interactions,  a  plane  free  shear  flow 
and  an  axisymmetric  jet  -  with  emphasis  on  acoustic  wave  propagation.  Our  numerical 
results  indicate  that  the  use  of  absorbing  layers  effectively  minimizes  numerical  reflection  in 
all  three  problems  considered. 
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1.  Introduction 

The  proper  treatment  of  computational  boundaries  is  crucial  for  any  numerical  solution  to 
a  set  of  partial  differential  equations  which  governs  fluid  motion  or  wave  propagation  in  a 
medium.  Various  techniques  have  been  developed  to  minimize  the  reflection  of  out-going 
waves.  A  review  can  be  found  in  Givoli  (1991).  Numerical  boundary  conditions  based  on 
the  characteristics  of  the  relevant  linearized  equations  and  their  asymptotic  solutions  in  the 
far  field  have  been  widely  used.  However,  such  boundary  conditions  are  not  satisfactory  if 
the  outflow  is  nonlinear  or  involves  multi-directional  waves.  As  a  possible  remedy,  a  buffer 
zone  abutting  the  computational  boundary,  in  which  the  governing  equations  are  modified, 
and  whose  role  is  to  absorb  the  incident  waves,  has  been  proposed.  In  this  buffer  zone,  the 
modifications  have  the  effect  of  either  removing  or  damping  reflected  waves  oriented  back 
towards  the  computational  domain.  Naturally,  the  buffer  zone  solutions  themselves  need  not 
necessarily  be  physical,  and  they  serve  only  to  prevent  contamination  of  the  solution  in  the 
physical  domain  of  interest  by  the  reflections  from  the  computational  boundaries.  Various 
types  of  buffer  zone  techniques  have  been  used  in  flow  simulations.  For  example,  Colonius 
et  al.(1993)  used  buffer  zones  in  which  the  solutions  were  filtered.  In  a  different  approach, 
Ta’asan  and  Nark  (1995)  modified  the  governing  equations  in  the  buffer  zone  to  change  the 
orientation  of  the  characteristics,  and  make  the  flow  supersonic  at  the  exit  plane.  Recently, 
Berenger(1994)  proposed  a  very  effective  Perfectly  Matched  Layer  technique  for  Maxwell’s 
equations.  In  this  approach,  the  equations  governing  the  so-called  matched  layer  are  split 
into  subcomponents  with  damping  terms  which  absorb  the  incident  waves  almost  perfectly. 
Following  Berenger,  Hu  (1996),  developed  an  analogous  technique  for  the  linearized  Euler 
equations,  and  provided  analytical  results  for  the  case  of  uniform  flow. 

In  this  paper,  we  follow  the  operator  splitting  principle  of  Berenger  (1994)  and  Hu  (1996) 
for  the  equations  governing  what  we  call  the  absorbing  layers  and  examine  their  effectiveness 
in  the  case  of  shock-vorticity  wave  interactions,  a  plane  free-shear  layer  and  an  axisymmetric 
jet.  The  emphasis  is  on  the  effectiveness  of  the  the  computational  boundaries  in  handling 
wave  propagation  including  sound  waves.  It  is  shown  that  the  absorbing  layer  technique  is 
very  effective  for  all  three  physical  problems.  The  next  section  describes  briefly  the  numerical 
models  used  in  this  study,  followed  by  the  section  on  results  and  conclusion. 


2.  Numerical  Models 

2.1  Shock  Wave  Interactions 

To  verify  the  applicability  of  the  absorbing  boundary  condition  technique  to  shock-turbulence 
and  shock-vortex  interaction  problems,  we  choose  the  numerical  model  of  Erlebacher,  Hus¬ 
sain!  and  Shu  (1997).  This  model  solves  the  fully  nonlinear  compressible  Euler  equations 
along  with  a  time  evolution  equation  for  the  shock  motion  for  the  purpose  of  fitting  the 
shock.  The  outflow  boundary  conditions  which  minimizes  wave  reflection  back  into  the  do¬ 
main  of  computation  are  of  crucial  importance  for  such  problems  as  they  involve  long-time 
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integrations.  The  present  case  focuses  on  the  interaction  of  a  single  vorticity  wave  with 
a  shock  wave,  and  the  results  of  course  carry  over  simply  to  a  randomly  distributed  wave 
system.  The  two  dimensional  Euler  equations  are  written  as 


dp  dp  dp  .du  dv. 


dv  dv  dv  dp 

dt  ^  dx  ^  dy  dy 

dp  ^  dp  _  dp  ,du  dv. 

The  computational  domain  has  the  shock  as  a  boundary  on  the  left  and  an  outflow 
boundary  on  the  right,  and  is  periodic  in  the  other  direction.  Fourth  order  Runge-Kutta 
scheme  is  used  for  time  integration,  and  the  spatial  derivatives  are  discretized  by  a  compact 
6**  order  scheme. 

In  the  absorbing  layer  at  the  right  boundary,  the  Euler  equations  are  split  into  a  lo¬ 
cally  one-dimensional  set  with  artiflcial  damping  terms.  Consider  the  pressure  equation,  for 
example,  in  computational  space: 


dp  dwi  8w2 


After  operator  splitting  and  addition  of  damping  terms,  the  pressure  equation  becomes 


dp2  dp  dw2 

aT  = W  ^'’aF  ■ 

in  the  absorbing  layer.  Here,  wi  and  W2  are  velocity  components  in  x-  and  y-directions, 
ai  and  02  are  contravariant  velocity  components  (which  include  the  effect  of  grid  motion) 
in  computational  space,  and  p  =  pi  -t-  ^2-  Locally  one-dimensional  equations  for  the  other 
variables  are  constructed  in  a  similar  manner.  The  damping  factor  ax  is  zero  in  a  layer 
parallel  to  the  X  direction;  similarly  ay  is  zero  in  a  layer  parallel  to  the  Y  direction  (see 
Figure  1).  However,  in  the  corner  region  both  these  damping  factors  are  positive. 
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2.2  Free  Shear  Layer 


In  order  to  evaluate  the  performance  of  the  absorbing-layer  technique  in  the  case  of  inviscid 
instability  waves,  we  solve  the  linearized  Euler  equations  in  a  Cartesian  {x,y)  coordinate 
system.  We  study  the  evolution  of  a  Kelvin-Helmholtz  instability  wave  as  it  propagates 
downstream  and  impinges  on  the  absorbing  layers.  In  this  case,  the  x-momentum  equation 
reads 


du  -du  dU  1  dp  ^ 


where 


U  =  ^  [{Ui  +  1/2)  +  {Ui  —  U2)  tanh(y)].  Absorbing  layers  are  used  at  the  upper,  lower  and 
right  boundaries.  Again,  the  afore-mentioned  operator  splitting  in  the  absorbing  layer  leads 
to  two  x-momentum  equations: 


dui 

'W 


frdu  1  dp 

+  (XxUi  -f-  U  — — h 

dx  pdx 


=  0 


du2  dU 


0 


where  u  =  ui  U2.  All  other  equations  are  treated  similarly.  These  equations  are  solved 
by  a  low-dissipation  and  low-dispersion  Runge-Kutta  scheme  which  is  formally  fourth  order 
accurate  (Hu,  Hussaini  and  Manthey,  1996). 


For  the  nonlinear  case  one  uses  again  an  approximate  time  independent  mean  flow  to 
split  the  Euler  equations  in  the  absorbing  layer.  Thus  the  stream-wise  velocity  for  two 
dimensional  flows  is  decomposed  into  three  components: 


u  =  U  +  ui  +  U2 


where  U  is  the  mean  velocity  as  in  the  linear  case.  Then  the  x-momentum  equation  is  written 
as 

du  1 

—  -t-  -f  VUy  -1-  -Px  =  0. 
dt  p 


This  equation  is  then  split  into  two  equations  as 

dui  1 

h  O'xUi  —  Px 

dt  p 


UUx  -f  -f  UUx 
P 


du2 

-h  (TyU2  =  -VUy 

All  other  equations  in  the  absorbing  layer  are  similarly  derived. 
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2.3  Axisymmetric  Jet 

The  compressible  axisymmetric  Euler  equations  for  the  jet  in  the  weak  conservation  form 
are  :  Qt  +  Gr  =  S,  where,  in  the  linearized  case, 
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In  the  above  equations,  p,  p,  rUz,  m^,  E  denote  the  fluctuating  components  of  pressure, 
density,  axial  and  radial  momentum,  total  energy  and  H  is  the  mean  enthalpy.  These  equa¬ 
tions  have  been  linearized  around  the  mean  velocity  {Ur,  Uz)  represented  by  an  error  function 
that  fits  experimental  measurements.  The  interior  equations  are  simply  split  into 

Qt  Fz  —  0,  Qt  +  Gr  =  S 

and  they  are  modified  in  the  absorbing  layer  as 

Ql  +  Fz  =  -azQ^  +  -SS  Q]  a  =  -OrQ'^  +  5^ 

where  Q  =  F  and  S  =  S^.  (We  used  =  0  in  this  study). 

In  the  nonlinear  case,  the  vectors  Q,S,F,  and  G  are  defined  as  follows. 
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We  use  the  fourth-order  MacCormack  method  which  has  been  successfully  used  in  earlier 
studies  by  Hayder  et  al.  (1996)  to  solve  the  linearized  Euler  equations,  and  by  Hayder  et  al. 
(1993)  and  Mankbadi  et  al.  (1994)  to  solve  the  Navier  Stokes  equations.  The  equations  are 
linearized  before  splitting  to  obtain  the  equations  for  the  absorbing  layer.  Thus,  we  get 

{Q^  -  Ql)t  +  {F-  Fo).  =  -  Ql) 


{Q^  —  Ql)t  +  {G  —  Go)r  =  —CTriQ^  —  Qo)  +  {S  —  So) 


where  the  subscript  0  denotes  mean  quantities. 


3.  Results 

In  the  case  of  the  shock-vorticity  wave  interaction,  we  consider  specifically  the  following 
simple  wave 

u 

u  —  Ui  =  eUi2y2-^cos[kxX  -f  kyy  —  Uit) 

k 

r-k 

V  =  —eUi2y2-^cos{kxX  +  kyy  —  Uit) 

k 

p=p=T=l 

as  the  upstream  condition  ahead  of  the  shock.  Ui  is  the  upstream  mean  velocity  normal  to 
the  undisturbed  shock,  ky  =  k  sin9,  kx  =  k  cos6  (where  k  is  upstream  wavenumber),  and 
e  =  0.001  measures  the  intensity  of  the  wave.  Our  standard  interior  domain  is  7.4  units  long 
with  185  uniformly  spaced  grid  points.  We  used  16  points  on  the  coordinate  axis  parallel 
to  the  shock.  An  absorbing  layer  abuts  the  right  outflow  boundary.  We  introduce  damping 
gradually  in  order  to  minimize  any  reflections  due  to  the  discretization  in  the  absorbing  layer. 
Unless  otherwise  mentioned,  we  use  0  =  30°,  K  =  2,  and  25  grid  points  (=  1  unit  in  length) 
in  the  buffer  layer  for  our  computations.  A  snapshot  of  pressure  in  the  interior  domain  at 
t  =  20  is  presented  in  Figure  2,  which  shows  how  well  the  out-going  waves  are  absorbed  with 
little  reflection.  To  measure  the  contamination  due  to  reflection,  the  solutions  are  compared 
with  a  reference  solution  obtained  by  computing  the  flow  in  a  much  larger  domain  with  the 
same  spatial  and  temporal  resolution.  We  follow  this  methodology  for  all  problems  in  this 
study.  Figure  3  compares  axial  variation  in  pressure  for  two  different  size  buffers  against 
the  large  domain  solution  at  t  =  20.  Because  of  modifications  to  the  governing  equations, 
the  solution  in  the  buffer  layer  is  irrelevant.  The  solution  in  the  interior  domain  for  a  buffer 
with  25  points  is  visually  indistinguishable  from  the  larger  domain  solution.  In  Figure  4, 
we  show  the  rms  error  (E)  in  pressure  at  the  ordinate  four  grid  points  upstream  of  the 
interface  between  the  computational  domain  and  the  absorbing  layer  as  a  function  of  the 
layer  thickness  measured  in  the  number  of  equidistant  points.  The  error  E  is  defined  as 

„  100 

|p;.„|v  N 
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where  p’’  is  the  pressure  from  the  reference  solution,  \Pmax\  maximum  amplitude  and  N 
is  the  number  of  grid  points  in  the  y-direction  (  A^=16  in  the  current  context).  E  measures 
the  numerical  error  in  the  solution,  which  includes  both  direct  and  induced  errors  due  to  the 
interaction  of  residual  reflections  from  the  outflow  boundary  with  the  flow  and  the  shock. 
As  expected,  E  decreases  as  the  layer  width  is  increased.  In  Figures  5  and  6,  we  show  the 
dependence  of  numerical  errors  on  the  angle  of  incidence  (0)  and  the  wave  number  (k).  The 
buffer  layer  is  more  effective  at  lower  incidence  angle  and  wavenumbers,  although  we  notice 
some  cross-overs  in  our  numerical  experiments.  At  later  times,  a  fraction  of  the  reflections 
from  the  outflow  boundary  propagates  upstream.  These  waves  can  then  reflect  back  and 
forth,  and  cause  what  we  call  induced  errors.  These  sometimes  constitute  a  significant 
portion  of  the  errors  shown  in  Figures  4-6  at  later  times. 

The  results  for  the  free-shear  layer  are  obtained  for  upper  and  lower  stream  mean  veloc¬ 
ities,  normalized  by  the  speed  of  sound,  equal  to  Ui  =  0.6  and  U2  =  0.2  respectively.  The 
eigenfunctions  of  the  Kelvin-Helmholtz  instability  wave  given  by  the  linear  stability  theory 
are  forced  at  the  inflow,  with  a  maximum  amplitude  e  equal  to  0.01.  We  solve  the  linearized 
Euler  equations  and  the  solution  agrees  with  the  linear  theory  very  well  in  eigenfunction 
and  growth  rate  comparisons.  In  Figure  7,  snapshots  of  axial  velocity  (Fig  7a)  and  pressure 
(Fig  7b)  are  shown,  and  in  Figure  8  we  present  the  amount  of  reflection  as  a  function  of 
the  layer  thickness.  We  observe  that  for  10  points  in  the  absorbing  layer,  the  amount  of 
reflection  (measured  four  grid  points  away  from  the  buffer  layer  boundary)  is  less  than  .03% 
of  the  amplitude  of  the  reference  pressure  fluctuation  from  the  large  domain  solution.  We 
also  solve  nonlinear  Euler  equations  where  the  nonlinearity  in  the  flow  is  significant.  The 
inflow  excitation  amplitude  (e)  is  kept  at  0.01,  but  the  interior  domain  is  three  times  longer. 
All  other  flow  parameters  are  the  same  as  in  the  linear  case.  The  error  in  pressure  four  grid 
points  away  from  the  buffer  layer  in  shown  in  Figure  9.  We  needed  a  larger  buffer  layer  for 
the  nonlinear  flow  simulations.  At  time  equal  to  3000,  errors  with  30  and  50  points  in  the 
buffer  layer  were  3.5%  and  4%  respectively.  Intuitively  one  expects  that  a  buffer  layer  to  be 
more  effective  if  nonlinear  effects  are  smaller.  This  may  be  the  principal  reason  for  larger 
errors  in  Figure  9  compared  to  Figure  8.  The  effect  of  nonlinearity  is  also  shown  in  Figure 

10,  where  we  compare  errors  for  simulations  with  two  different  levels  of  excitation  e  with  a 
buffer  of  50  grid  points. 

Finally,  for  the  case  of  the  excited  axisymmetric  jet,  we  assume  the  mean  Mach  number 
to  be  0.6.  At  the  inflow,  we  extrapolated  one  characteristic  variable  corresponding  to  the 
outgoing  acoustic  wave  from  the  interior  and  computed  the  other  three  characteristic  vari¬ 
ables  at  time  t  using  [/?,  u,u,p]  =  ei7e(qe®‘^*),  where  q  =  [p,u,v,p]  is  the  eigenfunction  given 
by  the  linear  stability  theory,  e  =  10“^,  u  =  1.05.  A  snapshot  of  pressure  is  shown  in  Figure 

11.  The  rms  pressure  error  (E)  in  the  immediate  neighborhood  (four  points  away  from  the 
buffer  layer)  of  the  layer  interface  is  plotted  in  Figure  12  for  time  equal  to  upto  50.  This 
error  becomes  quasi-periodic  and  the  maximum  error  for  25  grid  points  in  the  absorbing 
layer  is  about  0.015%.  Our  results  for  the  nonlinear  Euler  equations  are  shown  in  Figure  13. 
The  domain  size  is  10  units  long  for  both  linearized  and  nonlinear  Euler  simulation  of  the 
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excited  jet.  The  physical  parameters  are  the  same  for  both  the  linearized  and  the  nonlinear 
Euler  equations  for  the  jet  calculations. 


4.  Conclusions 

In  conclusion,  we  find  the  performance  of  the  absorbing-layer  technique  in  the  cases  of 
three  physical  problems  (using  three  different  numerical  algorithms)  is  quite  satisfactory. 
This  technique  offers  a  viable  alternative  to  the  traditional  boundary  treatments  based  on 
the  linearized  characteristics  or  asymptotic  solutions  in  the  far  field,  and  also  other  types  of 
buffer  layers.  It  also  promises  to  be  accurate  and  inexpensive  for  aeroacoustic  computations. 
Further  studies  are  warranted  to  put  this  methodology  on  a  firm  footing. 
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Figure  5:  Angle  dependence 
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Figure  6:  Wave  number  dependence 
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Figure  9:  Nonlinear  free  shear  layer 


Figure  10:  Effect  of  excitation  level 


Figure  11:  A  snapshot  of  pressure  in  the  jet 
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Figure  12:  Axisymmetric  jet  (Linearized  Euler) 


